Obtención de los datos

El objetivo principal de este estudio es analizar el efecto en la expresión genética debido a la infiltración mediante métodos en tiroides.

– Enlace al repositorio de github: https://github.com/MJoseom/ADO_RNAseq

El análisis de RNAseq se realiza a partir de los datos ya divididos en dos archivos: targets.csv y counts.csv empleando las herramientas de Bioconductor. Para realizar la práctica, se han seguido diferentes tutoriales [1] [2] [3] [4].

• Archivo targets.csv

Se trata del archivo de datos targets.csv con 292 observaciones, en este caso corresponde a datos de expresión (RNA-seq) pertenecientes a un análisis del tiroides. Para cada observación, se han anotado 9 variables o características. En este caso, se compara tres tipos de infiltración: 14 muestras pertenecientes al grupo Extensive lymphoid infiltrates (ELI), 236 muestras del grupo Not infiltrated tissues (NIT) y 42 muestras del grupo Small focal infiltrates (SFI). En la base de datos tenemos 0 missing values y 0 valores nulos. En este caso, las variables Grupo_analisis, body_site, molecular_data_type, sex, Group son consideradas factores.

Las variables o características consideradas son: Experiment, SRA_Sample, Sample_Name, Grupo_analisis, body_site, molecular_data_type, sex, Group, ShortName.

• Archivo counts.csv

Se trata del archivo de datos counts.csv con 292 muestras para las que se han analizado la expresión de 56202 genes en tiroides anotados con el Gencode ID como índice de cada fila. En este caso, el nombre de las columnas de la matriz Count coincide con la variable Sample_Name de los datos del archivo targets.csv. En la base de datos tenemos 0 missing values y 0 valores nulos. En este caso, las variables 292 son consideradas numéricas.

Selección de un subset

Se realiza una selección de 30 muestras de forma aleatoria, 10 para cada grupo de estudio (NIT, SFI, ELI).

1. Determinar la distribución de los grupos de estudio en función del sexo

2. Seleccionar 10 muestras de cada grupo de forma aleatoria del archivo targets.csv

set.seed(params$seed)
# Selección del subset
setTarget <- dataTarget %>% group_by(Group) %>% sample_n(10)
# Modificación la notación de la columna de muestras para evitar errores con la matriz de Count
setTarget$Sample_Name <- gsub("-", ".", setTarget$Sample_Name)
# Exportación la tabla con el subset de Targets
write.csv(setTarget, file= file.path(dataDir, "setTargets.csv"))

3. Analizar las características del dataSet

El dataset de Targets creado seleccionando muestras aleatorias de cada uno de los grupos, presenta 30 muestras para las que se han anotado 9 características o variables. En este caso, disponemos de tres tipos de infiltración: 10 muestras pertenecientes al grupo Extensive lymphoid infiltrates (ELI), 10 muestras del grupo Not infiltrated tissues (NIT) y 10 muestras del grupo Small focal infiltrates (SFI). En la base de datos tenemos 0 missing values y 0 valores nulos.

En este caso, dado que se seleccionan el mismo número de muestras para todos los grupos, no se mantien la misma relación de muestras en función del sexo y del grupo.

El identificador de las muestras seleccionadas para el dataset son:

4. Seleccionar los datos de expresión génica correspondientes al dataset del archivo counts.csv

Buscamos la coincidencia de los identificadores de la tabla 3 en las columnas las columnas correspondientes a las muestras del dataSet.

# Selección de las columnas con el identificador del dataset
setCount <- dataCount[,  sampleName]
# Eliminación de la versión en el nombre de cada gen con la notación de Ensambl
row.names(setCount) <- sapply(strsplit(row.names(setCount), split = ".", fixed=TRUE), function(a) a[1])
# Exportación la tabla con el subset de Counts
write.csv(setCount, file= file.path(dataDir, "setCounts.csv"))

El dataset de Counts creado presenta 30 muestras para las que se han analizado la expresión de 56202 genes en tiroides anotados con el Gencode ID como índice de cada fila. En este caso, el nombre de las columnas coincide con los 10 muestras de cada grupo (i.e., NIT, SFI, ELI) seleccionados de forma aleatoria en la Tabla 3. En la base de datos tenemos 0 missing values y 0 valores nulos.

Filtraje no específico

En la matriz setCount cada fila representa un gen con el código Ensembl y cada columna una librería de RNA secuenciada. Los valores corresponden a los datos sin procesar de fragmentos que se asignaron de forma única al gen respectivo en cada una de las librerías. Se empleará el paquete Deseq2 para realizar el análisis de expresión diferencial de los datos obtenidos por ultrasecuenciación. En este caso se determina como condición el efecto del tipo de infiltración (i.e., NIT, SFI, ELI) para cada una de las muestras.

dds <- DESeqDataSetFromMatrix(countData = setCount,
                                  colData = setTarget,
                                  design = ~ Group)
dds

## class: DESeqDataSet 
## dim: 56202 30 
## metadata(1): version
## assays(1): counts
## rownames(56202): ENSG00000223972 ENSG00000227232 ... ENSG00000210195
##   ENSG00000210196
## rowData names(0):
## colnames(30): GTEX.ZYY3.1926.SM.5GZXS GTEX.11XUK.0226.SM.5EQLW ...
##   GTEX.SIU8.0626.SM.2XCDN GTEX.TKQ1.0126.SM.33HB3
## colData names(9): Experiment SRA_Sample ... Group ShortName

En primer lugar, se realiza un filtro con el fin de eliminar todas aquellas filas que presentan únicamente 0 o bien un único valor de expresión genética. De este modo, se disminuye el tamaño de la matriz y aumentaremos la agilidad computacional. Así, pasamos de tener 56202 genes a reducir la cantidad de genes representativos hasta 43525.

Normalización

Con el objetivo de estabilizar la varianza entre los diferentes contajes de expresión génica, se ha aplicado el método VST (variance stabilizing transformation) y el rlog (regularized-logarithm transformation) del propio paquete DESseq2.

• VST

# VSD method
vsd <- vst(dds, blind = FALSE)
head(assay(vsd), 3)[,1:2]

##                 GTEX.ZYY3.1926.SM.5GZXS GTEX.11XUK.0226.SM.5EQLW
## ENSG00000223972                4.957860                 4.041418
## ENSG00000227232               10.206114                 8.957012
## ENSG00000243485                4.420783                 4.041418

head(colData(vsd), 3)

## DataFrame with 3 rows and 10 columns
##                           Experiment  SRA_Sample              Sample_Name
##                          <character> <character>              <character>
## GTEX.ZYY3.1926.SM.5GZXS    SRX568364   SRS627095  GTEX.ZYY3.1926.SM.5GZXS
## GTEX.11XUK.0226.SM.5EQLW   SRX619829   SRS644736 GTEX.11XUK.0226.SM.5EQLW
## GTEX.R55G.0726.SM.2TC6J    SRX204036   SRS374975  GTEX.R55G.0726.SM.2TC6J
##                          Grupo_analisis body_site        molecular_data_type
##                                <factor>  <factor>                   <factor>
## GTEX.ZYY3.1926.SM.5GZXS               3   Thyroid Allele-Specific Expression
## GTEX.11XUK.0226.SM.5EQLW              3   Thyroid              RNA Seq (NGS)
## GTEX.R55G.0726.SM.2TC6J               3   Thyroid              RNA Seq (NGS)
##                               sex    Group   ShortName        sizeFactor
##                          <factor> <factor> <character>         <numeric>
## GTEX.ZYY3.1926.SM.5GZXS    female      ELI   ZYY3-_ELI 0.873267401835296
## GTEX.11XUK.0226.SM.5EQLW   female      ELI   11XUK_ELI 0.901884946019608
## GTEX.R55G.0726.SM.2TC6J    female      ELI   R55G-_ELI  0.29590023913757

• rlog

# rlog method
rld <- rlog(dds, blind = FALSE)

## rlog() may take a few minutes with 30 or more samples,
## vst() is a much faster transformation

head(assay(rld), 3)[,1:2]

##                 GTEX.ZYY3.1926.SM.5GZXS GTEX.11XUK.0226.SM.5EQLW
## ENSG00000223972               1.8457820                1.0262641
## ENSG00000227232              10.0237124                9.0651040
## ENSG00000243485               0.3170474                0.1270095

Figure 1: Diagrama de dispersión de recuentos transformados de dos muestras por los métodos VST y rlog

La gráfica de dispersión de la Figura 1 nos muestra la transformación normalizada de counts en las dos primeras muestras empleando el método VST y rlog. Observamos que el rango de valores de rlog para el eje de las x es de -2.3 a 24.8, en cambio, el rango de valores para VST es menor para los valores inferiores, situándose en un rango de 4 a 21.7.

Un paso interesante en el análisis de RNAseq es determinar la similitud entre las muestras a partir de determinar la distancia entre éstas.

Para ello, realizamos una matriz de correlación de las distancias y con el fin de visualizarlas de una forma gráfica, realizamos un heatmap empleando el paquete pheatmap.

Figure 2: Heatmap de distancias de muestra a muestra usando los valores transformados de VST

En el heatmap de la Figura 2 observamos a partir del dendograma que hay 4 grandes clusters de muestras que presentan una cierta similitud. La distancia entre las muestras con mayor similitud, se encuentran con valores de 100 aproximadamente, tal y como se muestra en la paleta de colores del heatmap. Y éstas principalmente se concentran en los dos primeros clusters (primera rama izquierda del dendograma) pertenecientes a los grupos NIT y SFI. En general, las muestras de los diferentes grupos, se encuentran bastante agrupadas en el caso de la infiltración tipo ELI y la NIT. En cambio, las muestras del grupo SFI se encuentran más dispersas, mostrando mayor distancia entre ellas y por lo tanto, mayor variabilidad.

Figure 3: Visualización de las dos primeras componentes principales del PCA usando los valores transformados de VST

En el gráfico de la Figura 3 se han representado las dos primeras componentes del PCA (Principal Component Analysis) que explicarían el 72% de la variabilidad de las muestras, mayoritariamente por la PC1 que abarca el 64%. Observamos como en el caso anterior, que los grupos NIT y ELI, en menor medida, se encuentran agrupados. La variabilidad de la primera componente tiene una contribución muy alta del grupo NIT, que como podemos observar sitúa todas las muestras a la derecha de la gráfica, mientras que el grupo ELI se encuentra mayoritariamente a la izquierda. En cambio, las muestras de SFI se encuentran más repartidas por el centro del gráfico.

Figure 4: Visualización de las distancias con una matriz basada en MDS usando los valores transformados de VST

En el gráfico de la Figura 4, se muestra la distribución de los grupos empleando el Classical multidimensional scaling (MDS) para la matriz de distancias de los datos transformados de VST. El gráfico presenta un aspecto muy parecido al plot de PCA anterior (Figura 3), aunque en éste, además, se ha añadido la diferenciación por sexo, pero observamos que no hay un patrón claro entre las muestras de hombre o mujer.

Identificación de genes diferencialmente expresados

El parámetro log2FoldChange estima el cambio de expresión del gen debido al tratamiento, en este caso a los distintos métodos de infiltración. El parámetro padj corresponde al p-valor ajustado por el método Benjamini-Hochberg (BH). Con el fin de determinar que genes son signifcativos, podemos determinar un threshold para estos parámetros.

Con el fin de ser lo restrictivos, establecemos un threshold de 0.05 para el p-valor ajustado, mientras que dejamos el threshold de log2FoldChange por defecto.

Podemos realizar tres contrastes que comparen 2 a 2, los 3 grupos de infiltración en tiroides.

a) SFI-NIT

summary(SFIvsNIT)

## 
## out of 43520 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 588, 1.4%
## LFC < 0 (down)     : 72, 0.17%
## outliers [1]       : 0, 0%
## low counts [2]     : 16880, 39%
## (mean count < 3)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

En este caso observamos que los porcentajes de log2FoldChange para el resumen del análisis son bastante bajos, de modo que esperamos que la diferencia de tratamiento tenga poca incidencia en la expresión génica. Como hemos visto en la distribución de las Figuras 3 y 4, las muestras de SFI presentan mayor variabilidad y en consecuencia, se reduce el impacto en la expresión génica entre las dos condiciones comparadas. En este caso, aplicando los criterios establecidos, obtenemos 660 genes diferencialmente expresados y se descartarían 25985 genes.

En las Tablas 4 y 5 se destacan los 5 genes que en condiciones de infiltración SFI o NIT presentan un expresión más significativa, bien sea upregulated o downregulated. El contraste de significancia presentan un p-valor ajustado menor a 0.05 y se encuentran ordenadas de acuerdo al log2FoldChange. La primera columna de estas tablas indica el código del gen en GenCode.

b) ELI-NIT

summary(ELIvsNIT)

## 
## out of 43520 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 3394, 7.8%
## LFC < 0 (down)     : 1312, 3%
## outliers [1]       : 0, 0%
## low counts [2]     : 12661, 29%
## (mean count < 1)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

En este caso observamos que los porcentajes de log2FoldChange para el resumen del análisis son algo mayor que en en caso anterior, especialmente en la expresión de genes upregulated, de modo, que esperamos mayor número de genes significativos. En este caso, en las Figuras 3 y 4 veíamos que las muestras de estos grupos se encontraban en posiciones opuestas, de modo que podemos intuir que el método de infiltración NIT y ELI provocan mayores cambios en la expresión génica. En consecuencia, aplicando los criterios establecidos, obtenemos 4706 genes diferencialmente expresados y se descartarían 26158 genes.

En las Tablas 6 y 7 se destacan los 5 genes que en condiciones de infiltración ELI o NIT presentan un expresión más significativa, bien sea upregulated o downregulated. El contraste de significancia presentan un p-valor ajustado menor a 0.05 y se encuentran ordenadas de acuerdo al log2FoldChange. La primera columna de estas tablas indica el código del gen en GenCode.

c) ELI-SFI

summary(ELIvsSFI)

## 
## out of 43520 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 1324, 3%
## LFC < 0 (down)     : 295, 0.68%
## outliers [1]       : 0, 0%
## low counts [2]     : 14349, 33%
## (mean count < 2)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

En este caso observamos que los porcentajes de log2FoldChange para el resumen del análisis son igualmente bajos como en la primera comparación, de modo que esperamos que la diferencia de tratamiento tenga poca incidencia en la expresión génica. Como hemos visto en la distribución de las Figuras 3 y 4, las muestras de SFI presentan mayor variabilidad y en consecuencia, se reduce el impacto en la expresión génica entre las dos condiciones comparadas. En este caso, aplicando los criterios establecidos, obtenemos 1619 genes diferencialmente expresados y se descartarían 27557 genes.

En las Tablas 8 y 9 se destacan los 5 genes que en condiciones de infiltración ELI o SFI presentan un expresión más significativa, bien sea upregulated o downregulated. El contraste de significancia presentan un p-valor ajustado menor a 0.05 y se encuentran ordenadas de acuerdo al log2FoldChange. La primera columna de estas tablas indica el código del gen en GenCode.

Anotación de los resultados

• Plot del gen más significativo entre tratamientos

Figure 5: plotCounts del gen maś significativo en cada una de las comparaciones de los tratamientos

Figure 6: plotCounts del gen maś significativo en cada una de las comparaciones de los tratamientos

Figure 7: plotCounts del gen maś significativo en cada una de las comparaciones de los tratamientos

En las gráficas de plotCounts (Figuras 5, 6 y 7) se visibilizan los valores de counts para el gen más significativo, es decir, para el gen con el p-valor ajustado menor en la comparación de los grupos. En los 3 casos, destacan dos aspectos: a) los valores de ELI son más altos que en los otros dos grupos y los de NIT tienden a ser los más bajos y b) las mayores diferencias de valores se observan entre ELI y NIT, dado que SFI suele tener valores intermedios.

• Gene Clustering

Figure 8: Cluster de los genes con mayor variabilidad entre las muestras

En la Figura 8 se muestra la expresión de los 20 genes con mayor variabilidad entre las muestras que se encuentran anotadas con el grupo al que pertenecen. Se observa que una gran mayoría de genes, 14 de los 20, presentan una expresion downregulated especialmente en muestas del grupo NIT, mientras que para el resto se encuentran generalmente upregulated. En cambio, el primer cluster de genes (5 primeras filas) se expresan de forma diferencial en parte de las muestras NIT y SFI.

Análisis de significación biológicas

Una vez obtenemos la lista de genes expresados diferencialmente entre dos condiciones, debemos interpretar su relevancia biológica, es decir, conocer en que rutas metabólicas están implicados para conocer su función. Este análisis lo realizamos con la ayuda del paquete clusterProfiler que nos permite conocer dichas rutas de acuerdo al Entrez ID o al Symbol de cada gen.

En las siguientes Tablas 10, 11 y 12 se han anotado el Symboly el Entrez ID para cada uno de los genes significativamente diferenciados entre las distintas comparaciones. Aunque lamentablemente, algunos códigos no se encuentran en la base de datos de Bioconductor consultada (org.Hs.eg.db y EnsDb.Hsapiens.v86) probablemente dado que son transcritos.

a) SFI-NIT

Figure 9: Barplot para los términos de EnrichmentGO para el análisis del tipo de infiltración SFI vs NIT en la expresión génica

Figure 10: Red para los términos de EnrichmentGO para el análisis del tipo de infiltración SFI vs NIT en la expresión génica

En las Figura 9 y 10, observamos que gran parte de los genes diferenciados se encuentran involucrados en la unión a antígeno y al receptor de la immunoglobulina.

b) ELI-NIT

Figure 11: Barplot para los términos de EnrichmentGO para el análisis del tipo de infiltración ELI vs NIT en la expresión génica

Figure 12: Red para los términos de EnrichmentGO para el análisis del tipo de infiltración ELI vs NIT en la expresión génica

En las Figura 11 y 12, observamos que gran parte de los genes diferenciados se encuentran involucrados en la unión a antígeno y al receptor de la immunoglobulina. Además en este caso, también destacan con p-valores más altos pero significativos, los genes implicados en cell adhesion molecule binding, small GTPase binding y Ras GTPase binding.

c) ELI-SFI

Figure 13: Barplot para los términos de EnrichmentGO para el análisis del tipo de infiltración ELI vs SFI en la expresión génica

Figure 14: Red para los términos de EnrichmentGO para el análisis del tipo de infiltración ELI vs SFI en la expresión génica

En las Figura 13 y 14, observamos que en este caso, la función de los genes significativos se encuentra más repartida, destacando en este caso a diferencia de las comparaciones anteriores, nucleoside-triphosphatase regulator activity, DNA-binding transcription activator activity, RNA polymerase II-specific y protein serine/threonine kinase activity.

Abstract

El método de infiltración en tiroides se ha observado que puede afectar a la expresión genética. Para ello, se ha hecho un contaje de los tránscritos de 292 muestras, y se ha realizado un análisis de un subset de 30 de éstas distribuidas en los tres métodos de infiltración (i.e., SFI, NIT, ELI). Tras el análisis de la expresión de genes, se ha obtenido que hay una transcripción significativa principalmente entre los grupos NIT y ELI. Además, se ha observado que los genes se encuentra implicados en rutas metabólicas de unión a antígeno y al receptor de la immunoglobulina.

Introducción y Objetivo

El objetivo principal de este estudio es analizar el efecto en la expresión genética debido a la infiltración mediante métodos en tiroides.

No obstante, para una mejor contextualización de los resultados, sería necesario mayor información del estudio en el que se han recogido los datos analizados.

Materiales y Métodos

• Diseño del estudio

El estudio presenta un factor con 3 niveles, correspondiente a 3 métodos de infiltración en tiroides: a) Extensive lymphoid infiltrates (ELI); b) Not infiltrated tissues (NIT); c) Small focal infiltrates (SFI). La muestra original de participantes del estudio se ha reducido a una muestra de 30 sujetos, distribuyendo 10 réplicas para cada uno de los grupos. Además se han anotado otras características como el sexo, el experimento o un código identificativo.

Para cada una de las muestras, se ha realizado un contaje de la expresión transcripcional de una serie de genes implicados en tiroides, en total se han analizado 56202 genes.

• Diseño computacional

EL análisis de los datos se ha llevado a cabo empleando la versión 3.6.1 de R, así como funciones y paquetes pertenecientes al proyecto Bioconductor destinados a un análisis RNAseq. Con el fin de optimizar el protocolo, se han seguido los tutoriales presentados recientemente para dicho análisis [1] [2] [3] [4].

Resultados y Discusión

El análisis de los datos se ha llevado a cabo en primer lugar analizando la calidad de los datos. De modo que el primer paso es realizar un filtraje no específico eliminando aquellos genes para los que no se han obtenido contajes. De este filto, obtenemos 43525 genes que potencialmente pueden ser significativos. A continuación, se realiza una normalización de los datos usando el método variance stabilizing transformation (VST).

Seguidamente, se ha realizado un contraste en la que se ha especificado que los grupos a comparar son los tres métodos de infiltración, comparados dos a dos (SFI-NIT, ELI-NIT y ELI-SFI). Como resultado del test de significancia y aplicando cutoff para el p-valor ajustado con el método Benjamini and Hochberg menor a 0.05, se obtienen una serie de genes significativos para cada una de las comparaciones. Por último, se han analizado los términos de enrichment de GO database para agrupar los genes significativos en función de la ruta metabólica en la se ven implicada dichos genes.

De acuerdo con los resultados, se observa que hay mayores diferencias en la expresión de los grupos ELI y NIT, en cambio, las muestras de SFI presentan mayor variabilidad. Por lo tanto, se han obtenido 660, 4706 y 1619 genes significativos para las comparaciones SFI-NIT, ELI-NIT y ELI-SFI, respectivamente. La mayor parte de los genes significativos se encuentran involucrados en funciones muy diversas que incluyen la unión a antígeno y al receptor de la immunoglobulina, cell adhesion molecule binding, small GTPase binding, Ras GTPase binding, nucleoside-triphosphatase regulator activity o DNA-binding transcription activator activity.

Conclusión

Tras este informe, se ha conseguido analizar el contaje de genes transcritos aplicando un protocolo de RNAseq con el fin de valorar la respuesta génica en muestras con diferentes métodos de infiltración en tiroides. Se ha comprendido la importancia de cada paso del análisis y se han obtenido múltiples figuras que nos dan una imagen de la calidad de los datos, de la significancia de los test de comparación y de la participación de genes sigficativos en rutas metabólicas. No obstante, se requiere mayor información del estudio con el fin de comprender el contexto de éste. Así como también se requiere una mejor comprensión de toda la información extraída de cada gráfico y un análisis más detallado de cada uno de los genes significativos.